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ABSTRACT 

We present the simulation of 3D time dependent flow of rotating ideal gas falling 
into a Schwarzschild black hole. It is shown that also in the 3D case steady shocks 
are formed in a wide range of parameters (initial angular momentum and thermal 
energy) . Wc therefore highlight the stability of the phenomenon of shock formation in 
sub keplerian flows onto black holes, and reenforce the role of the shocks in the high 
luminosity emission from black hole candidates. The simulations have been performed 
using a parallelized code based on the Smoothed Particles Hydrodynamics method 
(SPH). We also discuss some properties of the shock problem that allow its use as a 
quantitative test of the accuracy of the used numerical method. This shows that the 
accuracy of SPH is acceptable although not excellent. 

Key words: accretion, accretion disks — black hole physics — hydrodynamics - 
instabilities 



1 INTRODUCTION 

As Bondi showed in his historical paper jBondilll952l) the 
occurrence of shocks in accretion phenomena on compact 
stellar objects is quite common. Since the radial speed can 
become supersonic the halting of the flow at a rigid wall 
leads necessarily to shock formation. 

Nevertheless the possibility of radial shock formation 
around black holes is less obvious since this compact object 
has no rigid wall. It is true that pre-heating of the infalling 
gas by radiation horn below or other non adiabatic pro- 
cesses can influence the properties of the falling gas and al- 
low shock solutions even in the spherically symmetric regime 
jBabul et al.lll989T) . 

However, if the gas is falling with some rotation, the 
centrifugal force can act as a "wall" and therefore trigger 
the shock formation even in the simple case of ideal gas. 
This fact was indeed suggested by the excellent pioneer gen- 
eral relativistic simu lations of Hawley, Smarr and Wilson 
jHawlev et al.lll984lh The question of the stability of the 
shock is critically related to the initial angular momentum 
of the gas, since for large values a standard keplerian disk 
will be formed, while for low values favorable conditions for 
standing shocks will occur. 

That such shocks can effectively stay at a fixed posi- 
tion has been p ut forward by S. Ch akrabarti (IChakrabartil 
1990), see also Chakr abarti jl998l) . It has been succes- 
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sively confirmed by the analytical work of Nakayama 
iNakavamalll994T) and by numeric al simulation experiments 
by C hakrabarti and Molteni llChakrabarti and Molteni 
1993) and following papers. The stability of the shock lo- 
cation is very relevant since, if it is stable, then the shock 
"engine" is always acting: it is not a transient phenomenon. 
Further studies, based on numerical simulations in 2D, axis 
symmetric configurations, revealed different kinds of oscil- 
lating behavio r over t he basic stable regime, s e e the results in 
IMolteni et all (Il996lh iMolteni et alJ (Il996bl) . iMolteni et alJ 
il999h . 

Recently the shock scenario has been shown to be very 
promising also to explain quasi perio dic oscillation proces ses 
(QPO) occurring in BH candidates iMolteni et alj|200ll) . 

Since all previous studies have been performed in 2D 
configurations (axial symmetric or in the Z=0 plane), it is 
natural to investigate the stability of the shocks in full 3D 
cases. Furthermore we add that this physical problem can 
be also considered a very critical test in the field of compu- 
tational fluid dynamics. 

The work proceeds as follows: in section 2 we resume 
the basic physics leading to shock formation and derive a 
simple formula to obtain the physical parameters necessary 
to generate the shocks, in section 3 we revise the numerical 
method and show how a test of any numerical code can be 
set up and in section 4 we show and discuss the new 3D 
simulation results. 
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2 SHOCK FORMATION IN UNVISCOUS 
SUB-KEPLERIAN FLOWS 

Let us resume the very basic physical ingredients of the phe- 
nomenon. We suppose that an unviscous gas is falling from 
very large distance onto a black hole. We will adopt the 
Paczyhski & Wiita potential jPaczvnski and Wiitalll980l) to 
mimic the Schwarzschild Black Hole force. It is given by : 

T / X G M * ,1\ 

*(»•) = (1) 

r-r g 

where r g — 2 ' G 2 M is the Schwarzschild radius. The basic 
physical effect can be easily understood in terms of "clas- 
sical" physics and it is well known that the Paczyriski & 
Wiita force reproduces many relativistic effects with high 
accuracy. 

To obtain the steady state solution we might integrate 
the differential equations for mass, momentum and energy 
conservations. In this case one has to start the space in- 
tegration fro m the sonic point as explained in the work by 
Chakrabarti (Chakrabarti 1990). However, in the case of un- 
viscous flow, it is easy to find an algebraic implicit solution, 
since the total energy is conserved (the Bernoulli theorem) 
and it can be fairly exploited to close the system of equations 
and find the solutions. 



2.1 ID Case 

Let us treat first the ID case. This derivation already ap- 
peared in the appendix of lMolteni et al.l il999T> . We repeat 
here the derivation pointing also to a peculiar aspect of the 
solution especially relevant to test code accuracy. By ID we 
mean that we are in an axially symmetric situation and the 
disk has no vertical extension [v z = 0, Z = 0). Note that 
the solution is obviously valid also in 2D dimensions with 
Z = (i.e. in the plane XY) if the symmetric condition is 
maintained. Both in the case of ID simulations in cylindrical 
coordinates and in 2D simulations (but with Z=0) the dif- 
ference between the theoretical shock position and the one 
resulting from simulations depends only on the numerical 
quality of the code; indeed there is no physical cause inter- 
vening to modify the results. Instead, in the 2D or 3D cases 
(with real vertical extension of the flow), the theoretical pro- 
cedure to derive the shock position needs the assumption of 
vertical equilibrium, that, as we will see, is not always sat- 
isfied. 

We assume also that the accreting gas is ideal (no vis- 
cosity and no cooling process is occurring). The angular mo- 
mentum per unit mass of the flow, A, must be conserved. In 
steady state regime we have also conservation of mass given 
by the equation: 



M 



-rpv r 



const 



(2) 



Here p is the surface density. In a conservative body 
force field with a potential VE' (r), also the pressure plays the 
role of a 'potential' and we have 

a(r) 2 



-v(rf+e(r)+ 



£M + ^ (r)= i„ (r)2+ __ 



+*(r) = 5(3) 



Here 7 is the adiabatic gas constant, B is the Bernoulli 
constant (the thermal energy at infinity), a is the sound 

speed a = \ xIL^ e j s the internal energy per unit mass, v 



is the speed of the flow and A is the angular momentum, so 
that \v 2 = \v\ + \vl = £ + \vl. 

If the entropy is conserved (no viscosity, no cooling, the 
presence of the shock will be taken into account few lines 
below) the polytropic relations are valid, P/p 1 = Po/Pq, 
so the relation between the density and sound speed p = 
K • a 2/(7 ~ 1) , where K = p /a 2 /h ~ 1) , is also valid. 

We have three unknown quantities p, v r , a, and three 
equations. The solutions are functions of r with angular mo- 
mentum and Bernoulli constant as parameters determining 
the specific shape of the solution. 

Using the radial Mach number m = — resolving a 
from the Bernoulli relation, and putting all terms in the 
continuity equation, we have the following implicit solution 
for the Mach number: 

M = r ■ m ■ a ■ K ■ oc f(m) ■ A (r, B, A) 

where / is function only of the mach number: 

f(m) = 



i 1 



1 1 2( 7 -l) 
7-1 J 



and A is function only of r (B and A are parameters) 



A (r, B, A) = r[B-V(r,X)]^=^ 



where 



GM* 



is the effective body force potential (gravitational plus cen- 
trifugal) . 

For a given M this equation can be easily solved. The 
/(m) function has a single maximum at m = 1, and it can 
be inverted both in the subsonic ^ m < 1 and supersonic 
m > 1 regions. The A function has in general two local 
minima at ri and V2 with n < r^, which can be determined 
by solving the algebraic equation dA/dr = numerically 
(high A values increase the centrifugal barrier and produce 
minima close to the BH). 

At large distances from the BH, the Mach number is 
m < 1, while close to the horizon it is m > 1. Since the 
/(m) • A (r, B, A) product must be constant along the flow, 
the maximum of f(m) must be at one of the minima of 
A(r), i.e. at n or 7-2. Therefore we can have two isentropic 
solutions mi(r) and mj(r) 



mi 



, 2 (r) = .r 1 



/(l)A A , s (n, 2 ) 



-4a 



(>•) 



(4) 



A standing shock can occur in the solution at r shock if 
the values m2(r s hock) > 1 and mi(r s (,«t) < 1 are related 
by the Hugoniot relation 



m 1 (r shock ) = h(m 2 {r shock )) 



(•>) 



where the Hugoniot relation for an hypothetical shock 
at a generic radius is given by 



h(m,2(r)) = 



2 + ( 7 -l)m 2 (r) 2 
277712 (r) 2 — (7 — 1) 



1/2 



(6) 



In general there can be two possible shock posi- 
tions, but only the outer one is stable as shown by 
the numerical simulations of Chakrabarti and Molteni 
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Figure 1. Analytical solution 



l|Chakrabarti and Moltenlll993l) and by the analytical work 
of Nakayama (Nakavama 1994). 

Figure 1 clarifies the procedure for a generic case. The 
supersonic line corresponds to the mi(r) solution while 
TO2(r) corresponds to the subsonic branch. 

2.2 Fictious and true rotation 

Let us note that the solution depends on the A function that 
contains the term: 



B 



A^ 
2r 2 



GM* 



(7) 



So, if we give to the flow an angular momentum lower 
than the theoretical one, but we add an equivalent centrifu- 
gal force, the resulting solution does not change. For exam- 
ple, we can obtain the same shock location if we give no 
real rotation to the flow and add to the gravitational force a 

,2 

fictitious force derived from the centrifugal potential , or 
viceversa, we can inject matter with an extremely fast real 
rotation, (even super keplerian!) and, at the same time, add 
a centripetal force. 

In this way, playing with real rotation and fictious cen- 
trifugal forces we can obtain a variety of solutions with the 
same radial shock position, same radial Mach number, same 
density profile. The important point is that the algebraic 
sum of real and fictious forces must be equal to the force 
required by the theory. It means that one of the following 
equations must be satisfied, according to the case 



theor 
y theor 



\eal + ^fict 



if 

if 



Xreal *C Xtheor 
Xreal Xtheor 



where Xtheor is the angular momentum value required 
by the theory to produce a shock at some specified radius. 
For consistency, if we define the real rotation as a fractio / 
of the theoretical one, i.e. X rea i = f -Xtheor, then the fictious 
force to be added is given by: 



Xfict 
Xfict 



Xtheor 
Xtheor 



VT 2 



p 



i 



if 
if 



/< i 
/>i 



It is clear that in the fast rotational cases the real shear 
motion can modify the results due to the numerical viscosity 
and other numerical code properties: accuracy, conservative 
structure, etc... 



We suggest the possibility to quantify the quality, Q, of 
a code as the ratio of the true rotation angular momentum 
to the theoretical one for which the shock is stable and close 
to the theoretical position: 

q = 1+ Xreal_ 
Atheor 

The obvious minimum requirement is that a code has 
Q=l so that, for no rotation, it produces the shock at the 
correct position. 



2.3 2D axis symmetric case 

For the 2D case (with real Z extension of the disk) the pro- 
cedure is similar. We maintain the hypothesis of axial sym- 
metry. We assume that the gas falls down in a condition 
of vertical equilibrium. Using the same variable names, the 
mass conservation equation is given by: 



M 



-4nrHpv r = const 



(8) 



where H is the half thickness of the disk. 

The vertical equilibrium hypothesis gives the well 
known expression for the half disk thickness: 



H = 



VGAL r (r - r g )a 



GM* 

The energy equation is still given by: 
A 2 1 2 , , P(r) 



(9) 



p{r) 



4-(r) 



B 



So again we have three unknown quantities p, v r , a, and 
three equations and we may find the solution exactly in the 
same way of the ID case. The difference will be in the exact 
values of the Mach function and of the geometrical function, 
A(r, X, B), but the shape are very similar: two minima in 
the geometrical function and one maximum in the Mach 
function. 

The Mach function is now given by: 
f ( m ) = — : — 



2 T 7-1J 



and the new A is function only of r (B and A are parameters) 



A (r, B, X) = r 5 ( r - r g ) [B — V (r, A)] <^-U 
where as previously stated 

w u A " GM * 

is the effective body force potential (gravitational plus cen- 
trifugal) . 



3 THE PHYSICAL MODEL AND THE 
NUMERICAL SIMULATIONS 

We performed 2D and 3D time dependent simulations of 
the motion of an ideal gas around a Schwarzschild black 
hole. The equations we integrate are in the lagrangian 
formulation (^ is the comoving derivative): 
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the continuity equation 

g = -pV-v (10) 

the momentum equation 

dv 1 

di p 

the energy equation 

df = P^ 
dt p 

the equation of state 

P = (7 - l)pe (13) 

Applying the SPH formalism fsee iMonaehanl jl992T0 . these 
are transformed into a set of ordinary differential equations 
(ODEs). These ODEs are subsequently integrated numeri- 
cally using standard methods. We use the usual summation 
form for the continuity equation: 



IV-v 



(11) 



(12) 



Pi = ^ m k Wik 



(14) 



(15) 



the momentum equation is given by: 
A N 

dt ^ Pi Pi 
the energy equation is: 

JV 

d ^ = -y mh (P^ + Pi+n ik )v ih -ViW ik (16) 

k — l 



Where Hit is the standard form of the artificial viscosity. 



3.1 Simulations parameters 

Here we show the results of 3D time dependent simulations 
and add also the results of 2D (Z=0) simulations to clarify 
the item of the code accuracy testing. We focus our atten- 
tion to a limited set of parameters which predicts a shock 
position close to the minimum radius allowed (around 5R g ) 
<Das et alJhoOll) . Therefore we are presenting the results 
concerning two critical cases. However we verified that also 
for other parameter values the shocks are formed and are 
stable, in good agreement with the predicted theoretical po- 
sition. The physical parameters of the cases we discuss here 
are given in Table 1 and they are well inside the stability 
region analyzed bv lDas et al.l i200ll) . 

The variables appearing in the Table 1 have the follow- 
ing meaning: ri is the inflow radius, v r i and a r i are the radial 
speed and sound speed at inflow radius, B is the Bernoulli 
constant, A is the angular momentum per unit mass, H is 
the vertical disk thickness at the inflow radius, K 3 hock is the 
predicted shock position, / is the fractio of the real rotation 
of the gas. 

The spatial domain (in cylindrical coordinates) is given 
by r = ~ r iy if — -r- 2tt , z — — 2r; -j- +2rj. 

The results have been obtained with a constant value of 
the spatial resolution h. This choice gives better results than 
those ones obtained with variable h; we think that the space 
varying h option has a larger amount of numerical viscosity 



due to the strong jump of h at the shock location. For 3D 
runs we used h = 0.5 and 7 = 4/3. 

The gas particles are injected (with the physical properties 
given in Table 1) from a set of injection points equally spaced 
onto the surface of a cylinder of radius n and height H, which 
is calculated from the vertical equilibrium assumption. A 
new particle is injected every time a sphere of radius h/2 
around each injector is found to be empty. 



3.2 Code parallelization 

Since our problem doesn't contain self gravitation, we used 
a spatial domain decomposition to parallelize the SPH code. 
We show here that this approach is very convenient for the 
kind of problem we are studying. Our parallel system con- 
sists of a biprocessor workstations cluster. As our proces- 
sors pool has local memory, we use the well known stan- 
dard: "Message Passing Interface" (MPI) to exchange data. 
Our physical problem (accretion disks around central star) 
has a space domain with some degree of central symmetry, 
so it seems logical to divide the overall domain into sub 
domains having this type of symmetry. The parallelization 
paradigm we use is the Multiple Instruction Multiple Data 
one (MIMD). With this decomposition and paradigms, the 
code performing the calculations for one sub domain is the 
standard serial code used for one single domain, with minor 
modifications. Information exchanges occur always between 
two consecutive domains. 

The three dimensional computational domain is decom- 
posed into concentrically cylindrical coronae (sub domains) 
with centers coincident with the coordinate origin. Figure 2 
shows an horizontal cut view, in a case in which the domain 
is decomposed into four sub domains. All the SPH particles, 
having vector radius projection (in XY plane) lying between 
two consecutive cylinders, are assigned to the same sub do- 
main. The number of sub domains is equal or multiple of 
the number of processors in the cluster. 

As the SPH particles have a short-range interaction, 
only the particles lying within a distance 2h from the 
limiting cylinders interact with their analogue particles lay- 
ing in the previous cylinder or in the next one. 

In any moment each processor manages m internal par- 
ticles plus nbi n edge particles from the inner edge region and 
nb ou t edge particles from the outer edge region. 

In each iteration loop the m particle of a sub domain are 
evolved taking into account the presence of the all nb par- 
ticles (for density, pressure forces etc. calculation). Instead 
the nb particles are evolved according to the kinematic and 
dynamic properties they had as internal particles of their 
original domain. At the end of the time step and before to 
start the next one a check is performed to verify if some 
internal particle has migrated into another domain. In this 
case the particles data are exchanged. 



4 RESULTS AND DATA ANALYSIS 

We show the basic results of the simulations using one panel 
for each case. In each panel we have 4 figures: the figure top- 
left shows the number of particles versus time; the figure 
top-right shows the radial Mach number of all the parti- 
cles folded, in cp, around the z axis; the figure bottom-left 
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Panel 1. Simulation with f=0.0 
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Figure 2. domain subdivision 



shows the radial speed; the figure bottom-right shows the 
distribution of the angular momentum of the particles, the 
scattering of the values is just due to the conservation of the 
total angular momentum. 

As it is clear from Figure 1.1 the system reaches the 
equilibrium. Figure 1.2, with the folded radial Mach number 
values, shows that the shock is located at Rshock = 10.5, 
a position close to the theoretical one (Rshock = 6.3). A 
not perfect agreement between the theoretical shock position 
and the simulation result is not surprising: a discrepancy 
has to be expected due to the vertical motion of the gas in 
3D geometry, while the theory assumes vertical equilibrium 
(the Mach number profile is also very close to the theoretical 
one). Figure 1.4 shows that also the mean value of A is close 
to the real zero value. 

We then simulated the same flow but increasing the 
amount of true rotation. For a fractio / = 0.1 the flow still 
reaches the equilibrium state, but the shock is located some- 
what outwards compared to the theoretical position. 

For / = 0.2 the results are shown in Panel 3. The shock 
is steady, but at a larger radius. 

For / = 0.4, see Panel 4, the system is still in steady 
state and the shock is further shifted outwards. It is clearly 
seen that the number of particles versus time is increased. 

For / ^ 0.7 (not shown) the shock position is unstable: 
it travels steadily outwards. The number of particles is al- 
ways increasing. We stopped the simulation before the shock 
reached the outer boundary of our integration domain. 

In general the figures showing the radial speed versus 
r point out that, as the true rotation increases, the number 
of outflowing particles increases. We inject with a fixed 
angular momentum, but the numerical code viscosity, 
while conserving the total angular momentum, produces 
an angular momentum redistribution and therefore the 
particles with large angular momentum flow away. 
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Figure 1.1 Number of particles versus time 




Figure 1.2 Radial Mach Number versus R 




Figure 1.3 Radial speed versus R 




Figure 1.4 True angular momentum versus R 
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Panel 2. Simulation with f=0.1 
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Figure 2.1 Number of particles versus time 




Figure 2.2 Radial Mach Number versus R 



We made also 2D simulations, on the Z = plane. In 
this case the agreement between the simulation results and 
the theory should be perfect; the discrepancies are fully at- 
tributable to the code properties. We checked that the 2D 
SPH simulations with different f values show clearly the ef- 
fect of the numerical shear viscosity. Panel 5 shows the re- 
sults obtained with / = 0.5. The shock is stable. In figure 5.1 
the shock position versus time is plotted, in figure 5.2 the ra- 
dial Mach number (folded in <p) is plotted versus R, in figure 

5.3 the radial speed (folded in y>) is displayed and in in figure 

5.4 the particle distribution in the XY plane is presented. 
In this case the agreement between theory and simulation 
is quite good. Note that the symmetry of the flow is so well 
conserved that the folded data appear as one dimensional 
in R. The physical parameters are E — 0.002, A = 1.8. For 
gas with 7 = 5/3 the shock must be Rshock = 6.78, in good 
agreement with the simulation results. The oscillations in 
Figure 5.1 are due to the algorithm of the shock identifica- 
tion. The SPH numerical parameters are h = 0.2 and we use 
a = 1, ft = 2 for the artificial viscosity coefficients. In the 
case with f=0.75, here not shown, the shock is still stable, 
but the position is different from the predicted one. In gen- 
eral, if we compare 2D and 3D simulations with the same 
parameters and resolution, the 3D case is more stable than 
the 2D one. We interpret this fact as due to the impossibil- 
ity of ejecting (along Z) high angular momentum gas outside 
the domain in the 2D case, while in the 3D case the high 
angular momentum particles are vertically ejected leading 
to a stabilizing effect. 
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Figure 2.3 Radial speed versus R 




Figure 2.4 True angular momentum versus R 



5 CONCLUSIONS 

We showed that even in 3D time dependent simulations 
the standing shocks predicted by Chakrabarti's analysis are 
formed. In the numerical simulations, for values of the pa- 
rameters in the theoretical stable zone, the shocks are stable. 
This fact confirms the relevance of the shock phenomena in 
accretion flows around black holes in galactic and extra- 
galactic sources. Indeed the standard keplerian disk models 
around galactic black hole candidates suffer instabilities that 
make not completely clear the physical scenario. In the case 
of the extragalactic sources like AGN (active galactic nu- 
clei), supposed to contain black holes, while the keplerian 
disk can easily take account of the low energy emission, the 
high luminosity at high energies, X and 7, is still an open 
problem. 

Due to the high CPU cost of 3D simulations the vari- 
ability properties of the shock flow exhibited in the reduced 
dimensional simulations (2D in r-z axis symmetric configu- 
ration and 2D XY plane configuration) will be the subject 
of a subsequent paper. 

We showed also that this physical problem can be used 
as a test for any fluid dynamics code since it is possible 
to derive the analytical profile of the solution of the flow 
and compare it with the simulated one for different degrees 
of real rotation of the fluid. In this way it is possible to 
measure the numerical shear viscosity of the code. We show 
that the 3D SPH code detects the shock at a good level of 
accuracy. 
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Panel 3. Simulation with f=0.2 Panel 4. Simulation with f=0.4 
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Figure 3.1 Number of particles versus time Figure 4.1 Number of particles versus time 




Figure 3.2 Radial Mach Number versus R Figure 4.2 Radial Mach Number versus R 




Figure 3.3 Radial speed versus R Figure 4.3 Radial speed versus R 




Figure 3.4 True angular momentum versus R 



Figure 4.4 True angular momentum versus R 
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Panel 5. Simulation with f=0.4 
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Figure 5.1 Shock position versus time 
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Figure 5.2 Radial Mach Number versus R 
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Figure 5.3 Radial speed versus R 



10 15 20 



Figure 5.4 Particle distribution in the XY plane 
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Table 1. 





r; 






A 


H 


f 


shock 


B 


I 


38 


-0.076 


0.069 


1.625 


16 


0.0 


5.76 


0.0045 


II 


50 


-0.072 


0.058 


1.65 


20 


0.1 


6.32 


0.003 


III 


50 


-0.072 


0.058 


1.65 


20 


0.2 


6.32 


0.003 


IV 


50 


-0.072 


0.058 


1.65 


20 


0.4 


6.32 


0.003 



